METODO DE RUNGE-KUTTA DE CUARTO ORDEN PARA LA RESOLUCION NUMERICA DEL MODELO NEURONAL DE FITZHUGH-NAGUNO.

Clear[f, x, y, desde, hasta, n, h, i, nodo, inicio1, inicio2, L1, L2, 
      L3, L4, k1, k2, k3, k4]
f[x_, y_, z_] := 0.23 - y*(y - 0.2)*(y - 1) - z
g[x_, y_, z_] := 0.02*(y - 0.5*z)
desde := 0. ;
hasta := 300
inicio1 = {0};
inicio2 = {0};
n = 300;
h = ((hasta - desde)/n);
nodo = Table[desde + i h, {i, 0, n}];
For[i = 2, 
        i ≤ n + 1, i++, 
k1 = 
          f[nodo[[i - 1]], 
            inicio1[[i - 1]], 
            inicio2[[i - 
                1]]]; 
L1 = 
          g[nodo[[i - 1]], 
            inicio1[[i - 1]], 
            inicio2[[i - 
                1]]]; 
k2 = 
          f[nodo[[i - 1]] + (h/2), 
            inicio1[[i - 
                  1]] + ((h k1)/2), 
            inicio2[[i - 
                  1]] + ((h L1)/2)]; 
L2 = 
          g[nodo[[i - 1]] + (h/2), 
            inicio1[[i - 
                  1]] + ((h k1)/2), 
            inicio2[[i - 
                  1]] + ((h L1)/2)]; 
k3 = 
          f[nodo[[i - 1]] + (h/2), 
            inicio1[[i - 
                  1]] + ((h k2)/2), 
            inicio2[[i - 
                  1]] + ((h L2)/2)]; 
L3 = 
          g[nodo[[i - 1]] + (h/2), 
            inicio1[[i - 
                  1]] + ((h k2)/2), 
            inicio2[[i - 
                  1]] + ((h L2)/2)]; 
k4 = 
          f[nodo[[i - 1]] + h, 
            inicio1[[i - 1]] + h k3, 
            inicio2[[i - 1]] + 
              h L3]; 
L4 = 
          g[nodo[[i - 1]] + h, 
            inicio1[[i - 1]] + h k3, 
            inicio2[[i - 1]] + 
              h L3]; 
AppendTo[inicio1, 
          inicio1[[i - 1]] + 
            (1/6) h (k1 + 2 k2 + 2 k3 + 
                  k4)]; 
AppendTo[inicio2, 
          inicio2[[i - 1]] + 
            (1/6) h (L1 + 2 L2 + 2 L3 + L4)];];

Representación gráfica de la solución aproximada correspondiente al potencial de membrana

graf1 = 
        ListPlot[
          Table[{nodo[[i]], 
              inicio1[[i]]}, {i, 
              n + 1}], PlotStyle -> RGBColor[1, 0, 0]];

Representación gráfica de la solución aproximada correspondiente a l conductancia de los iones

graf2 = 
        ListPlot[
          Table[{nodo[[i]], 
              inicio2[[i]]}, {i, 
              n + 1}], 
          PlotStyle -> 
            RGBColor[0, 0, 1]];

Superposición de ambas soluciones

Show[graf1, graf2]

Plano fase

ListPlot[
        Table[{inicio1[[i]], 
            inicio2[[i]]}, {i, n + 1}], 
        PlotStyle -> RGBColor[0, 0, 1], 
        AspectRatio -> 1 ];

El comportamiento a "largo plazo" es el siguiente: las órbitas giran en torno al punto de equilibrio (0.11 , 022). Es decir, las dos variables se comportan de forma periódica.